home *** CD-ROM | disk | FTP | other *** search
/ The Best of Down Under Games / The Best of Down Under Games.iso / 3dfx Screen Savers / VoodooLights / Sources / mat.c < prev    next >
C/C++ Source or Header  |  1997-07-12  |  10KB  |  390 lines
  1. /*------------------------------------------------------/
  2. /                                                        /
  3. /    Copyright 1997, SΘrgio Durte <smd@di.fct.unl.pt>    /
  4. /                                                        /
  5. /------------------------------------------------------*/
  6.  
  7. #include <stdio.h>
  8. #include <math.h>
  9. #include <stdlib.h>
  10.  
  11. #include "defines.h"
  12. #include "mat.h"
  13.  
  14. XYZ xyz_Zero = { 0.0, 0.0, 0.0 } ;
  15.  
  16. XYZW xyzw_Zero = { 0.0, 0.0, 0.0, 1.0 } ;
  17.  
  18. int mat_ipow( int b,  int e )
  19. {
  20.   int result = 1 ;
  21.   while( e-- ) result *= b ;
  22.   
  23.   return result ;
  24. }
  25.  
  26. void mat_scalev( Float s, XYZ *in, XYZ *out)
  27. { 
  28.   register Float f = s ;
  29.   out->x = in->x * f ;
  30.   out->y = in->y * f ;
  31.   out->z = in->z * f ;
  32. }
  33.  
  34. Float mat_dotp( XYZ *a, XYZ *b)
  35. {
  36.   return ((a->x) * (b->x) + (a->y) * (b->y) + (a->z) * (b->z)) ;
  37. }
  38.  
  39. void mat_addv( XYZ *a, XYZ *b, XYZ *out )
  40. { 
  41.   out->x =a->x + b->x ;
  42.   out->y =a->y + b->y ;
  43.   out->z =a->z + b->z ;
  44. }
  45.  
  46. void mat_subv( XYZ *a, XYZ *b, XYZ *out )
  47. { 
  48.   out->x =a->x - b->x ;
  49.   out->y =a->y - b->y ;
  50.   out->z =a->z - b->z ;
  51. }
  52.  
  53. void mat_combv( Float t, XYZ *a, XYZ *b, XYZ *out )
  54. {
  55.   out->x = t * a->x + b->x ;
  56.   out->y = t * a->y + b->y ;
  57.   out->z = t * a->z + b->z ;
  58. }
  59.  
  60. Float mat_length( XYZ *a )
  61. { 
  62.   return  (Float) sqrt ( sqr(a->x ) + sqr(a->y) + sqr(a->z) ) ;
  63. }
  64.  
  65. Float mat_normalize( XYZ *a )
  66. {
  67.    Float n, nn ;
  68.  
  69.    nn = sqr(a->x) + sqr(a->y) + sqr(a->z) ; 
  70.    if ( nn < ZERO ) return 0.0 ; 
  71.  
  72.    nn = (Float) sqrt(nn) ;
  73.    n = 1.0f / nn ;
  74.    a->x *= n ; 
  75.    a->y *= n ; 
  76.    a->z *= n ; 
  77.    return nn ;
  78. }   
  79. void mat_crossp( XYZ *a, XYZ *b, XYZ *out )
  80. { 
  81.   out->x = a->y * b->z - b->y * a->z ;
  82.   out->y = b->x * a->z - a->x * b->z ;
  83.   out->z = a->x * b->y - b->x * a->y ; 
  84. }
  85.  
  86.  
  87. void mat_rotv( XYZ *e, XYZ *p, XYZ *o )
  88. {
  89.   Float m = (Float)sqrt( sqr( p->y ) + sqr( p->z ) ) ;
  90.   Float im = 1.0f / m ;
  91.  
  92.   o->x = m * e->x + p->x * e->z ;
  93.   o->y = ( -p->x * p->y * e->x  + p->z * e->y ) * im + p->y * e->z ;
  94.   o->z = ( -p->x * p->z * e->x  - p->y * e->y ) * im + p->z * e->z ;
  95. }
  96.  
  97. void mat_randomDir( XYZ *d )
  98. {
  99.     d->x = 1 - 2 * rnd() ;
  100.     d->y = 1 - 2 * rnd() ;
  101.     d->z = 1 - 2 * rnd() ;
  102.     mat_direction( d ) ;
  103. }
  104.  
  105. void mat_printv( XYZ *a )
  106. {
  107.   fprintf(stderr,">%12g\t%12g\t%12g\n",(double)a->x , (double)a->y , (double)a->z ) ;
  108. }
  109.  
  110. void mat_printXYZW( XYZW *a )
  111. {
  112.     Float x, y, z ;
  113.     
  114.     x = (Float)(a->x / a->w) ;
  115.     y = (Float)(a->y / a->w) ;
  116.     z = (Float)(a->z / a->w) ;
  117.     
  118.     fprintf(stderr,">%8g\t%8g\t%8g\t%8g\n",(double)a->x , (double)a->y , (double)a->z, (double)a->w) ;
  119.     fprintf(stderr,"3D ->[%8g\t%8g\t%8g]\n", (double)x, (double)y, (double)z ) ;
  120. }
  121.  
  122. void mat_mult_m_XYZp( XYZ *v, Matrix a, XYZ *out )
  123. {
  124.   out->x = v->x * a[0][0] + v->y * a[0][1] + v->z * a[0][2] + a[0][3] ;
  125.   out->y = v->x * a[1][0] + v->y * a[1][1] + v->z * a[1][2] + a[1][3] ;
  126.   out->z = v->x * a[2][0] + v->y * a[2][1] + v->z * a[2][2] + a[2][3] ;
  127. }
  128.  
  129. void mat_mult_m_XYZWp( XYZ *v, Matrix a, XYZW *out )
  130. {
  131.   out->x = v->x * a[0][0] + v->y * a[0][1] + v->z * a[0][2] + a[0][3] ;
  132.   out->y = v->x * a[1][0] + v->y * a[1][1] + v->z * a[1][2] + a[1][3] ;
  133.   out->z = v->x * a[2][0] + v->y * a[2][1] + v->z * a[2][2] + a[2][3] ;
  134.   out->w = v->x * a[3][0] + v->y * a[3][1] + v->z * a[3][2] + a[3][3] ;
  135. }
  136.  
  137. void mat_mult_m_XYZv( XYZ *v, Matrix a, XYZ *out )
  138. {
  139.   out->x = v->x * a[0][0] + v->y * a[0][1] + v->z * a[0][2] ;
  140.   out->y = v->x * a[1][0] + v->y * a[1][1] + v->z * a[1][2] ;
  141.   out->z = v->x * a[2][0] + v->y * a[2][1] + v->z * a[2][2] ;
  142. }
  143.  
  144. void mat_multpm( XYZ *v, Matrix a, XYZW *out )
  145. {
  146.   out->x = v->x * a[0][0] + v->y * a[0][1] + v->z * a[0][2] + a[0][3] ;
  147.   out->y = v->x * a[1][0] + v->y * a[1][1] + v->z * a[1][2] + a[1][3] ;
  148.   out->z = v->x * a[2][0] + v->y * a[2][1] + v->z * a[2][2] + a[2][3] ;
  149.   out->w = v->x * a[3][0] + v->y * a[3][1] + v->z * a[3][2] + a[3][3] ;
  150. }
  151.  
  152. void mat_multvm(XYZ *v, Matrix a, XYZ *out )
  153. { 
  154.   out->x = v->x * a[0][0] + v->y * a[0][1] + v->z * a[0][2] ;
  155.   out->y = v->x * a[1][0] + v->y * a[1][1] + v->z * a[1][2] ;
  156.   out->z = v->x * a[2][0] + v->y * a[2][1] + v->z * a[2][2] ;
  157. }
  158.  
  159. void mat_mult_tm_XYZv(XYZ *v, Matrix a, XYZ *out )
  160. { 
  161.   out->x = v->x * a[0][0] + v->y * a[1][0] + v->z * a[2][0] ;
  162.   out->y = v->x * a[0][1] + v->y * a[1][1] + v->z * a[2][1] ;
  163.   out->z = v->x * a[0][2] + v->y * a[1][2] + v->z * a[2][2] ;
  164. }
  165. void mat_mult_tm_XYZp(XYZ *v, Matrix a, XYZ *out )
  166. {  
  167.   out->x = v->x * a[0][0] + v->y * a[1][0] + v->z * a[2][0] ;
  168.   out->y = v->x * a[0][1] + v->y * a[1][1] + v->z * a[2][1] ;
  169.   out->z = v->x * a[0][2] + v->y * a[1][2] + v->z * a[2][2] ;
  170.  
  171.   out->x -= a[0][0] * a[0][3] + a[1][0] * a[1][3] + a[2][0] * a[2][3] ;
  172.   out->y -= a[0][1] * a[0][3] + a[1][1] * a[1][3] + a[2][1] * a[2][3] ;
  173.   out->z -= a[0][2] * a[0][3] + a[1][2] * a[1][3] + a[2][2] * a[2][3] ;
  174. }
  175.  
  176. void mat_multvtm(XYZ *v, Matrix a, XYZ *out )
  177. { 
  178.   out->x = v->x * a[0][0] + v->y * a[1][0] + v->z * a[2][0] ;
  179.   out->y = v->x * a[0][1] + v->y * a[1][1] + v->z * a[2][1] ;
  180.   out->z = v->x * a[0][2] + v->y * a[1][2] + v->z * a[2][2] ;
  181. }
  182.  
  183. void mat_multpim(XYZ *v, Matrix a, XYZ *out )
  184. {  
  185.   out->x = v->x * a[0][0] + v->y * a[1][0] + v->z * a[2][0] ;
  186.   out->y = v->x * a[0][1] + v->y * a[1][1] + v->z * a[2][1] ;
  187.   out->z = v->x * a[0][2] + v->y * a[1][2] + v->z * a[2][2] ;
  188.  
  189.   out->x -= a[0][0] * a[0][3] + a[1][0] * a[1][3] + a[2][0] * a[2][3] ;
  190.   out->y -= a[0][1] * a[0][3] + a[1][1] * a[1][3] + a[2][1] * a[2][3] ;
  191.   out->z -= a[0][2] * a[0][3] + a[1][2] * a[1][3] + a[2][2] * a[2][3] ;
  192. }
  193.  
  194. void mat_summ( Matrix a, Matrix b, Matrix out )
  195. {
  196.   register int k;
  197.   for(k=3 ; k>=0 ; k--) 
  198.     {
  199.       out[0][k] = a[0][k] + b[0][k] ;
  200.       out[1][k] = a[1][k] + b[1][k] ;
  201.       out[2][k] = a[2][k] + b[2][k] ;
  202.       out[3][k] = a[3][k] + b[3][k] ;
  203.     }
  204. }
  205.  
  206. void mat_subm( Matrix a, Matrix b, Matrix out)
  207. {
  208.   register int k;
  209.   for(k=3 ; k>=0 ; k--) 
  210.     {
  211.       out[0][k] = a[0][k] - b[0][k] ;
  212.       out[1][k] = a[1][k] - b[1][k] ;
  213.       out[2][k] = a[2][k] - b[2][k] ;
  214.       out[3][k] = a[3][k] - b[3][k] ;
  215.     }
  216. }
  217.  
  218. void mat_multm(Matrix a, Matrix b, Matrix out )
  219. {
  220.   register int i, j ;
  221.   for(i=3 ; i>=0 ; i--)
  222.       for(j=3 ; j>=0 ; j--)
  223.         out[i][j] = a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j] + a[i][3] * b[3][j] ; 
  224. }
  225.  
  226. void mat_idm( Matrix out )
  227. { 
  228.  register int i, j ;
  229.  for( i = 0 ; i < 4 ; i++ )
  230.      for( j = 0 ; j < 4 ; j++ ) out[i][j] = ( i == j ? 1.0f : 0.0f ) ;
  231. }
  232.  
  233. void mat_transl( XYZ *v, Matrix out )
  234. {
  235.   mat_idm( out ) ;
  236.   out[0][3] = v->x ;
  237.   out[1][3] = v->y ;
  238.   out[2][3] = v->z ;
  239. }
  240.  
  241. void mat_scale( XYZ *v, Matrix out )
  242. {
  243.   mat_idm( out );
  244.   out[0][0] = v->x ;
  245.   out[1][1] = v->y ;
  246.   out[2][2] = v->z ;
  247. }
  248.  
  249. void mat_rotx(Float degrees, Matrix out )
  250. { 
  251.   Float c, s ;
  252.   c = (Float)cos( RAD( degrees ) ) ;
  253.   s = (Float)sin( RAD( degrees) ) ;
  254.   mat_idm( out ) ;
  255.   out[1][1] = c ;
  256.   out[2][2] = c ;
  257.   out[1][2] = -s ;
  258.   out[2][1] = s ;
  259. }
  260.  
  261. void mat_roty(Float degrees, Matrix out ) 
  262. { 
  263.   Float c, s ;
  264.   c = (Float)cos( RAD( degrees ) ) ;
  265.   s = (Float)sin( RAD( degrees ) ) ;
  266.   mat_idm( out ) ;
  267.   out[0][0] = c ;
  268.   out[2][2] = c ;
  269.   out[0][2] = s ;
  270.   out[2][0] = -s ;
  271. }
  272.  
  273. void mat_rotz(Float degrees, Matrix out ) 
  274. { 
  275.   Float c,s;
  276.   c = (Float)cos( RAD( degrees ) ) ;
  277.   s = (Float)sin( RAD( degrees ) ) ;
  278.   mat_idm( out ) ;
  279.   out[0][0] = c ;
  280.   out[1][1] = c ;
  281.   out[0][1] = -s ;
  282.   out[1][0] = s ;
  283. }
  284.  
  285. void mat_shearx( Float sy, Float sz, Matrix out )
  286. {
  287.     mat_idm( out ) ;
  288.     out[1][0] = sy ;
  289.     out[2][0] = sz ;
  290. }
  291. void mat_sheary( Float sx, Float sz, Matrix out )
  292. {
  293.     mat_idm( out ) ;
  294.     out[0][1] = sx ;
  295.     out[2][1] = sz ;
  296. }
  297.  
  298. void mat_shearz( Float sx, Float sy, Matrix out )
  299. {
  300.     mat_idm( out ) ;
  301.     out[0][2] = sx ;
  302.     out[1][2] = sy ;
  303. }
  304.  
  305.  
  306. void mat_transpose( Matrix in, Matrix out )
  307. {
  308.     int i, j ;
  309.     for( i = 0 ; i < 4 ; i++ )
  310.         for( j = 0 ; j < 4 ; j++ ) out[i][j] = in[j][i] ;
  311. }
  312.  
  313. void mat_invm( Matrix m, Matrix out )
  314. {
  315.   register int i, j ;
  316.  
  317.   mat_idm( out ) ;
  318.   for ( i = 0 ; i < 3 ; i++ )
  319.     for ( j = 0 ; j < 3; j++ )
  320.       out[i][j] = m[j][i] ;
  321.   
  322.   for ( i = 0 ; i < 3 ; i++ )
  323.     out[i][3] = -( m[0][i] * m[0][3] + m[1][i] * m[1][3] + m[2][i] * m[2][3] ) ;
  324.   
  325. }
  326.  
  327.  
  328. void mat_printm( char *s, Matrix m )
  329. {
  330.   int j ;
  331.   fprintf(stderr, "%s\n", s) ;
  332.   for( j = 0 ; j < 4 ; j++ )
  333.     fprintf(stderr, "%8.6g\t%8.6g\t%8.6g\t%8.6g\n", m[j][0], m[j][1], m[j][2], m[j][3] ) ;
  334. }
  335.  
  336.  
  337. #define mat_Accumulate    if( temp >= 0.0 ) pos += temp ; else neg += temp 
  338. #define Precision_Limit   1e-15 
  339.  
  340. Bool mat_ainvm( Matrix i , Matrix o )
  341. {
  342.   double det_1 ;
  343.   double pos, neg, temp ;
  344.  
  345.   pos = neg = 0.0 ;
  346.   temp = i[0][0] * i[1][1] * i[2][2] ;
  347.   mat_Accumulate ;
  348.   temp = i[0][1] * i[1][2] * i[2][0] ;
  349.   mat_Accumulate ;
  350.   temp = i[0][2] * i[1][0] * i[2][1] ;
  351.   mat_Accumulate ;
  352.   temp = -i[0][2] * i[1][1] * i[2][0] ;
  353.   mat_Accumulate ;
  354.   temp = -i[0][1] * i[1][0] * i[2][2] ;
  355.   mat_Accumulate ;
  356.   temp = -i[0][0] * i[1][2] * i[2][1] ;
  357.   mat_Accumulate ;
  358.   det_1 = pos + neg ;
  359.  
  360.   if( det_1 == 0.0 || fabs( det_1 / ( pos - neg )) < Precision_Limit )
  361.     {
  362.       fprintf( stderr," Singular Matrix, cannot invert.\n");
  363.       return False ;
  364.     }
  365.  
  366.   det_1 = 1.0 / det_1 ;
  367.   
  368.   o[0][0] = (Float)( ( i[1][1] * i[2][2] - i[1][2] * i[2][1] ) * det_1 ) ;
  369.   o[1][0] = (Float)(-( i[1][0] * i[2][2] - i[1][2] * i[2][0] ) * det_1 ) ;
  370.   o[2][0] = (Float)( ( i[1][0] * i[2][1] - i[1][1] * i[2][0] ) * det_1 ) ;
  371.   o[3][0] = 0.0 ;
  372.  
  373.   o[0][1] = (Float)(-( i[0][1] * i[2][2] - i[0][2] * i[2][1] ) * det_1 ) ;
  374.   o[1][1] = (Float)( ( i[0][0] * i[2][2] - i[0][2] * i[2][0] ) * det_1 ) ;
  375.   o[2][1] = (Float)(-( i[0][0] * i[2][1] - i[0][1] * i[2][0] ) * det_1 ) ;
  376.   o[3][1] = 0.0 ;
  377.  
  378.   o[0][2] = (Float)( ( i[0][1] * i[1][2] - i[0][2] * i[1][1] ) * det_1 ) ;
  379.   o[1][2] = (Float)(-( i[0][0] * i[1][2] - i[0][2] * i[1][0] ) * det_1 ) ;
  380.   o[2][2] = (Float)( ( i[0][0] * i[1][1] - i[0][1] * i[1][0] ) * det_1 ) ;  
  381.   o[3][2] = 0.0 ;
  382.  
  383.   o[0][3] = (Float)(- ( i[0][3] * o[0][0] + i[1][3] * o[0][1] + i[2][3] * o[0][2] )) ;
  384.   o[1][3] = (Float)(- ( i[0][3] * o[1][0] + i[1][3] * o[1][1] + i[2][3] * o[1][2] )) ;
  385.   o[2][3] = (Float)(- ( i[0][3] * o[2][0] + i[1][3] * o[2][1] + i[2][3] * o[2][2] )) ;
  386.   o[3][3] = 1.0 ;
  387.   
  388.   return True ;
  389. }
  390.